function [teststat,pval,yfit,pfit]=testjfe(ybar,phat,n,vary,covyd,vard,K,alpha,makegraph,maxchunksize)
% inputs:
%   ybar: J-vector of mean outcomes by judge
%   phat: J-vector of estimated propensities by judge
%   n: J-vector of number of cases by judge
%   vary: J-vector of variance of y by judge
%   covyd: J-vector of covariance between y and d by judge
%   vard: J-vector of variance of d by judge
%   K: absolute value bound on treatment effect magnitude
%   graph: indicator if you want to graph best constrained fit
%
% output:   1 = reject 
%           0 = fail to reject



% trim empirical phats away from zero or one
phat=max(phat,1-alpha.^(1./n));
phat=min(phat,alpha.^(1./n));
% Number of judges
J = length(ybar);

if nargin<10 maxchunksize=J;
end
maxchunksize=maxchunksize
% sort by phat
%origorder=(1:J)';
[phat,ind]=sort(phat);
ybar=ybar(ind);
n=n(ind);
vary=vary(ind);
covyd=covyd(ind);
vard=vard(ind);
%origorder=origorder(ind);

% set seed
rng(1);
% changing from random chunks to propensity-based chunks:
%[~,randorder]=sort(rand(J,1));
%chunks=ceil(randorder/maxchunksize)
chunks=ceil((1:J)'/maxchunksize);

numchunks=max(chunks);
teststats=zeros(numchunks,1);
pfitmax=0;
yfitmax=0;
phatmax=0;
ybarmax=0;
yfit=zeros(J,1);
pfit=yfit;


        
% impose constraints within chunks
for ch=1:numchunks
   %fprintf('on chunk %d of %d\n',ch,numchunks);
    ybarch=ybar(chunks==ch);
    Jch=length(ybarch);
    phatch=phat(chunks==ch);
    if Jch==1
        yfit(chunks==ch)=ybarch;
        pfit(chunks==ch)=phatch;
    else
        clear model;
        nch=n(chunks==ch);
        varych=vary(chunks==ch);
        covydch=covyd(chunks==ch);
        vardch=vard(chunks==ch);
        
        % I am organizing the choice variables as follows. We have 6 choice
        % variables per judge:
%         1:  yj
%         2:  pj
%         3:  T1j
%         4:  T2j
%         5:  abs(T1j)
%         6:  abs(T2j)
        % So there are 6*J of these
        % and we also have choice variables associated with the (j,j')
        % constraints:
%         1:  yj-yj'
%         2:  pj-pj'
%         3:  abs(yj-yj')
%         4:  abs(pj-pj')
%         There are 4*J*(J-1)/2=2*J*(J-1) of these
% And finally one element to be the max of the studentized deviations
       % So the choice variable vector has 6*J + 2*J*(J-1)+1 rows (assuming
       % we impose all of the pairwise constraints.)

        % The linear equality constraint matrix will be:
%         Aeq = [A1;
%              A2],
%          where A1 define quantities within each judge only, and A2 defines 
%          quantities associated with (j,j') constraints.
%         A1 and A2 each have 6*J + 2*J*(J-1)+1 columns (one for each choice variable).
%         A1 has 2*J (because we define T1j and T2j for each judge) rows and A2 has 2*J*(J-1)/2 = J*(J-1) rows.
        A1 = sparse(2*Jch,6*Jch+2*Jch*(Jch-1)+1);
        A2 = sparse(Jch*(Jch-1),6*Jch+2*Jch*(Jch-1)+1);

        % The linear inequality constraint matrix will have J*(J-1)/2 rows,
        % because we only have one for each pair of judges:
        Aineq=sparse(Jch*(Jch-1)/2,6*Jch+2*Jch*(Jch-1)+1);
        
        % The linear constraints need a rhs vector, which will have
        % 2*J+3/2*J*(J-1) rows.
        rhs = zeros(2*Jch+3/2*Jch*(Jch-1),1);
        
        % We need to collect the indices for abs(T1j) and abs(T2j):
        absts = sort([5:6:6*Jch-1,6:6:6*Jch]');
        % We need a vector for the objective function (all zeros except for a one in the last spot):
        f=[zeros(6*Jch+2*Jch*(Jch-1),1);1];

%         We also need absolute value constraints. We have 2*J  for defining
%         abs(T1j) and abs(T2j), and J*(J-1) for defining
%         abs(yj-yjp) and abs(pj-pjp).
        % Let's first build the constraints for each judge
        for j = 1:Jch
           rhoj=covydch(j)/sqrt(varych(j)*vardch(j));
           yjind = 6*j-5;
           pjind = 6*j-4;
           T1ind = 6*j-3;
           T2ind = 6*j-2;
           absT1ind = 6*j-1;
           absT2ind = 6*j;
           % define T1j = (ybarj - yj - K(phatj - pj))/se1, so
           % yj - K*pj + se1*T1j = ybarj - K*phatj
           se1=sqrt(varych(j)+K^2*vardch(j)-2*K*covydch(j))/sqrt(nch(j));
           A1(2*j-1,yjind) = 1;
           A1(2*j-1,pjind) = -K;
           A1(2*j-1,T1ind) = se1;
           rhs(2*j-1) = ybarch(j)-K*phatch(j);
           % define T2j = (ybarj - yj + K(phatj - pj))/se2, so
           % yj + K*pj + se2*T2j = ybarj + K*phatj
           se2=sqrt(varych(j)+K^2*vardch(j)+2*K*covydch(j))/sqrt(nch(j));
           A1(2*j,yjind) = 1;
           A1(2*j,pjind) = K;
           A1(2*j,T2ind) = se2;
           rhs(2*j) = ybarch(j)+K*phatch(j);
           % define abs(T1j)
           model.genconabs(2*j-1).argvar=T1ind;
           model.genconabs(2*j-1).resvar=absT1ind;
           % and abs(T2j)
           model.genconabs(2*j).argvar=T2ind;
           model.genconabs(2*j).resvar=absT2ind;  
        end
        
        % Now build the pairwise constraints
        for j = 1:Jch-1
            %fprintf('building constraints, j = %d of %d\n',j,Jch-1);
            for jp=j+1:Jch
                % define counter for the current pair
                cc=Jch*(j-1)-j*(j+1)/2+jp;
                % define choice vector index for this pair's quantities:
                yjind =6*j -5;
                yjpind=6*jp-5;
                pjind =6*j -4;
                pjpind=6*jp-4;
                ydiffind=6*Jch+4*cc-3;
                pdiffind=6*Jch+4*cc-2;
                absydiffind=6*Jch+4*cc-1;
                abspdiffind=6*Jch+4*cc;
                
                % Linear equality constraints
                % first yj - yjp
                % yj
                A2(2*cc-1,yjind)=1;
                % yjp
                A2(2*cc-1,yjpind)=-1;
                % yj - yjp
                A2(2*cc-1,ydiffind)=-1;
                
                % Now pj -pjp
                % pj
                A2(2*cc,pjind)=1;
                % pjp
                A2(2*cc,pjpind)=-1;
                % pj - pjp
                A2(2*cc,pdiffind)=-1;
                
                % Nonlinear (abs) constraints
                % abs(yj-yjp)
                model.genconabs(2*Jch+2*cc-1).argvar=ydiffind;
                model.genconabs(2*Jch+2*cc-1).resvar=absydiffind;
                % abs(pj-pjp)
                model.genconabs(2*Jch+2*cc).argvar=pdiffind;
                model.genconabs(2*Jch+2*cc).resvar=abspdiffind;
                    
                % Linear inequality (slope) constraints abs(a(j,jp))-K*abs(b(j,jp))<=0
                Aineq(cc,absydiffind)=1;
                Aineq(cc,abspdiffind)=-K;
            end
        end
        % Now define max of the studentized deviations
        model.genconmax(1).resvar = 6*Jch+2*Jch*(Jch-1)+1;
        model.genconmax(1).vars=absts;
        
        % Now put things together
        senseeq=repmat('=',2*Jch+Jch*(Jch-1),1);
        senseineq=repmat('<',Jch*(Jch-1)/2,1);
        sense=[senseeq;senseineq];
             
        Aeq=[A1;A2];
        A=[Aeq;Aineq];
        
        % now bounds:
        % yj in (-inf,inf)
        % pj in (0,1)
        % T1j in (-inf,inf)
        % T2j in (-inf,inf)
        % |T1j| in (0,inf)
        % |T2j| in (0,inf)
        
        % yj - yjp in (-inf,inf)
        % pj - pjp in (-1,1)
        % |yj - yjp| in (0,inf)
        % |pj - pjp| in (0,1)
        
        % max{Ts} in (0,inf)
        lbj = [-inf;0;-inf;-inf;0;0];
        ubj = [inf;1;inf;inf;inf;inf];
        lbjjp = [-inf;-1;0;0];
        ubjjp = [inf;1;inf;1];
        
        lb = [repmat(lbj,Jch,1);repmat(lbjjp,Jch*(Jch-1)/2,1);0];
        ub = [repmat(ubj,Jch,1);repmat(ubjjp,Jch*(Jch-1)/2,1);inf];
        % Set objective: x
        model.obj = f;
        model.A   = A;
        model.rhs=rhs;
        model.sense=sense;
        model.lb=lb;
        model.ub=ub;

        %fprintf('solving model\n');
        result = gurobi(model);
        teststats(ch)=result.objval;
        xch=result.x;
        yfitch=xch(1:6:6*Jch-5);
        pfitch=xch(2:6:6*Jch-4);       
        yfit(chunks==ch)=yfitch;
        pfit(chunks==ch)=pfitch;
    end
end
teststat=max(teststats);

pvala=1;
for j = 1:J
   rhoj=(vary(j)-K^2*vard(j))/sqrt((vary(j)+K^2*vard(j))^2-(2*K*covyd(j))^2);
   sj = [1,rhoj;rhoj,1];
   %Sigma(2*j-1:2*j,2*j-1:2*j)=sj;
   f1=mvncdf([teststat;teststat],[0;0],sj);
   f2=mvncdf([-teststat;-teststat],[0;0],sj);
   f3=mvncdf([-teststat;teststat],[0;0],sj);
   if isnan(f3) f3=0;end
   Fj = f1+f2-2*f3;
   pvala = pvala*Fj;
end

pval = 1-pvala;

if makegraph>0
    [pfit,ind]=sort(pfit);
    yfit=yfit(ind);
    phat=phat(ind);
    ybar=ybar(ind);
    %Plot results
    %testgraph=figure('Position',[520 378 1120 420],...
     %   'PaperOrientation','landscape',...
      %  'PaperPosition',[1,3,9,4]);
    testgraph = figure('PaperOrientation','portrait','PaperPosition',[1,4,6,6]);
    figure(testgraph);
    scatter(phat,ybar);
    title(sprintf('Average Outcomes and Propensity by Judge\n'));
    xlabel('propensity');
    ylabel('average outcome');
    minp=min([pfit;phat]);
    maxp=max([pfit;phat]);
    miny=min([yfit;ybar]);
    maxy=max([yfit;ybar]);
    diffp=maxp-minp;
    diffy=maxy-miny;
    maxdiff=max(diffp,diffy);
    diffdiffp=maxdiff-diffp;
    diffdiffy=maxdiff-diffy;
    %axis([minp minp+maxdiff miny miny+maxdiff]);
    axis([minp-diffdiffp/2 minp+diffp+diffdiffp/2 miny-diffdiffy/2 miny+diffy+diffdiffy/2]);
    hold on
    quiver(pfit,yfit,phat-pfit,ybar-yfit,0,'ShowArrowHead','off','Color',[1 .5 0])
    plot(pfit,yfit,'Color',[1 0 0],'LineWidth',2);
    pbaspect([1 1 1]);
    hold off
end
end